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Abstract 

We consider the problem of constructing Lyapunov functions for linear differential equations with delays. For such systems it 
is known that exponential stability implies the existence of a positive Lyapunov function which is quadratic on the space of 
continuous functions. We give an explicit parametrization of a sequence of finite-dimensional subsets of the cone of positive 
Lyapunov functions using positive semidefinite matrices. This allows stability analysis of linear time-delay systems to be 
formulated as a semidefinite program. 
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1 Introduction 

In this paper we present an approach to the construction of Lyapunov functions for systems with time-delays. 
Specifically, we are interested in systems of the form 



±{t) = ^2 A i x it ~ hi) 



where x(t) s K ra . In the simplest case we are given the delays ho, . . . ,hk and the matrices Aq, . . . , Ak and we would 
like to determine whether the system is stable. For such systems it is known that if the system is stable, then there 
exists a Lyapunov function of the form 



V(<f>) 
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c/>(s) T N(s,t)(t)(t) dsdt, 



where M and N are piecewise-continuous matrix-valued functions, and h = m&x{ho, . . . , hk}. Here (j> : [—h, 0] — > M n 
is an element of the state space, which in this case is the space of continuous functions mapping [— h, 0] to R n . 
The function V is thus a quadratic form on the state space. The derivative is also such a quadratic form, and the 
matrix- valued functions which define it depend affinely on M and N. 



* Corresponding author M. Peet. 

Email addresses: matthew.peet@inria.fr (Matthew M. Peet), antonis@eng.ox.ac.uk (Antonis Papachristodoulou), 
lall@stanford.edu (Sanjay Lall). 



Preprint submitted to Automatica 



1 February 2008 



In this paper we develop an approach which uses semidefinite programming to construct piecewise-continuous func- 
tions M and N such that the function V is positive and its derivative is negative. The functions we construct are 
piecewise-polynomial. Roughly speaking, we show that 



<f>(a) 



M(s) 



0(0) 
0(a) 



ds 



is positive for all (f> if and only if there exists a piecewise continuous matrix-valued function T such that 



M(t) 



T(t) 




> for all t 



[ T(t) dt = 0. 

J-h 



That is, we convert positivity of the integral to pointwise positivity of M. This result is stated precisely in Theorem 5. 
Pointwise positivity may then be easily enforced, and in the case of positivity on the real line this is equivalent to a 
sum-of-squares constraint. The constraint that T integrates to zero is a simple linear constraint on the coefficients 
of T. Notice that the sufficient condition that M(s) be pointwise nonnegative is conservative, and as the equivalence 
above shows it is easy to generate examples where V\ is nonnegative even though M(s) is not pointwise nonnegative. 



We also give a necessary and sufficient condition for positivity of 




0(s) T N(s,t) <j>{t) dsdt. 



Roughly speaking, if N is polynomial, then the given form is positive if and only if there exists a positive semidefinite 
matrix Q y such that N(s,t) — Z(s) T QZ(t), where Z is a vector of monomials. This condition is expressed as a 
constraint on the coefficients of N and is stated in Theorem 9. Note that pointwise positivity of N is not sufficient 
for positivity of the functional. The condition that the derivative of the Lyapunov function be negative is similarly 
enforced. 



1.1 Background 



The use of Lyapunov functions on an infinite dimensional space to analyze differential equations with delay originates 
with the work of Krasovskii 7]. For linear systems, quadratic Lyapunov functions were first considered by Repin 
[lH . The book of Gu, Kharitonov, and Chen m| presents many useful results in this area, and further references may 
be found there as well as in Hale and Lunel Q , Kolmanovskii and Myshkis @ and Niculescu @ . The idea of using 
sum-of-squares polynomials together with semidefinite programming to construct Lyapunov functions originates 
in Parrilo @. 



1.2 Notation 



Let N denote the set of nonnegative integers. Let S n be the set ofnxn real symmetric matrices, and for X € S" we 
write X y to mean that X is positive semidefinite. For two matrices A,B, we denote the Kronecker product by 
A® B. For X any Banach space and I C M any interval, let 0(7, X) be the space of all functions 

n(I,X) = {f:I^X} 

and let C(I, X) be the Banach space of bounded continuous functions 

C(I . X) ={/:/—> X | / is continuous and bounded } 

equipped with the norm 

ll/H =sup||/(t)||x. 
tei 
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We will omit the range space when it is clear from the context; for example we write C[a, b] to mean C([a,b], X). 
A function is called C n (I,X) if the i th derivative exists and is a continuous function for i = 0, . . . , n. A function 
/ G C[a, b] is called piecewise continuous if there exists a Unite number of points a < fti < ■ ■ ■ < hk < b such that / 
is continuous at all x G [a, b]\{h\, . . . , hk} and its right and left-hand limits exist at {fti, . . . , hk}. 

Define also the projection H t : Q[— h, oo) — > f2[— h, 0] for t > and h > by 

(H t x)(s) = x(t + s) for all s G [-h,0]. 

We follow the usual convention and denote H t x by x t . 



2 System Formulation 



Suppose = h < hi < ■ ■ ■ < hk — h. Define the sets H = {— fto, . . . , —hk} and H c = [—h,0]\H and suppose 
A ,...,4 ei" x ". We consider linear differential equations with delay, of the form 

k 

X (t) = ^2 A i X ( f - h i) for a11 1 > °! C 1 ) 

where the trajectory x : [—h, oo) — > R n . The boundary conditions are specified by a given function <f> : \—h, 0] — > M" 
and the constraint 

x(t) = 0(i) for alH e [-/i,0]. (2) 

If G C[— h, 0], then there exists a unique function a; satisfying (1) and (2). The system is called exponentially 
stable if there exists u > and ael such that for every initial condition <f> E C[— h, 0] the corresponding solution 
x satisfies 

||a;(t)|| < oe _<rt ||0|| for alii > 0. 

We write the solution as an explicit function of the initial conditions using the map G : C[— h,Q] — > Q[—h,oo), 
defined by 

(G<j))(t) = x(t) for alii > -ft, 

where a; is the unique solution of (1) and (2) corresponding to initial condition cf>. Also for s > define the /kw map 
r s : C[-h, 0] -» C[-ft, 0] by 

which maps the state of the system x t to the state at a later time x t + s — F s Xt- 



2.1 Lyapunov Functions 



Suppose V : C[— ft, 0] — > K. We use the notion of derivative as follows. Define the Lie derivative of V with respect 
to T by 

V{<j>) = lim sup - (V(r r <f>) - V(4>)) . 

We will use the notation V for both the Lie derivative and the usual derivative, and state explicitly which we mean if 
it is not clear from context. We will consider the set X of quadratic functions, where V G X if there exist piecewise 
continuous functions M : [-ft, 0) -> § 2n and N : [-ft, 0) x [-ft, 0) -> R nxn such that 







T 

M(s) 


0(0)" 


l-h 


>( S )_ 





,0 ,0 

rfs+ / / (f){s) T N{s,t)(j)(t)dsdt. 
J-h J-h 



(3) 



The following important result shows that for linear systems with delay, the system is exponentially stable if and 
only if there exists a quadratic Lyapunov function. 
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Theorem 1. The linear system defined by equations (1) and (2) is exponentially stable if and only if there exists a 
Lie-differentiable function V G X and e > such that for all <p G C[—h, 0] 



V(4>) >e||0(O)|| 2 
V(4>) < -eU(0)\\ 2 - 



(4) 



Further V G X may be chosen such that the corresponding functions M and N of equation (3) have the following 
smoothness property: M(s) and N(s,t) are bounded and continuous at all s,t such that s G H c and t G H°. 

Proof. Sec Kharitonov and Hinrichsen [4| for a recent proof. □ 

3 Positivity of Integrals 

The goal of this section to present results which enable us to computationally find functions V G X which satisfy 
the positivity conditions in (4) and have the form 



V(y) = f 

J-h 



M(t) 



y(o) 
.»(*). 



dt. 



Before stating the main result in Theorem 5, we give a few preliminary lemmas. 

Lemma 2. Suppose f : [—h,0] — ► M is piecewise continuous. Then the following are equivalent. 



(i) / f(t)dt>a 

J-h 

(ii) There exists a function g: [—h, 0] — > M which is piecewise continuous and satisfies 

f(t)+g(t)>0 for all t, 



g(t) dt = 0. 

Proof. The direction (ii) (i) is immediate. To show the other direction, suppose (i) holds, and let g be 

9(t) = -/(*) + \ f h f{s) ds for all t. 

Then g satisfies (ii). 



□ 



The second lemma shows that minimizing over continuous functions is as good as minimizing over piecewise contin- 
uous functions. 

Lemma 3. Suppose H = {— ho, . . . , — hk} and let H° = [—h,0]\H. Let f : [—h,0] x R n — > K be continuous on 
H c x W 1 , and suppose there exists a bounded function z : [— h, 0] — > R, continuous on H c , such that for all t G [—h, 0] 

f(t,z(t))=Mf(t,x) 
Further suppose for each bounded set X C K™ the set 

{f(t,x) | xex,t g [-h,o]} 

is bounded. Then 



/o M 
f(t,y(t))dt= / iid f(t,x)dt 
-h J-h x 



(5) 
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Proof. Let 

K= I mif{t,x)dt 
J-h x 

It is easy to see that 

inf / f(t,y(t))dt> K 
since if not there would exist some continuous function y and some interval on which 

f{t,y(t)) <m£f(t,x) 

which is clearly impossible. 

We now show that the left-hand side of (5) is also less than or equal to K, and hence equals K. We need to show 
that for any e > there exists y G C[— h, 0] such that 



[ f(t,y(t)) 
J-h 



dt < K + s. 



To do this, for each n € N define the set H n C R by 

fe-i 



H n = [J {hi - a/n, hi + a/n) 



i=i 



and choose a > sufficiently small so that Hi C (—h, 0). Let z be as in the hypothesis of the lemma, and pick M 
and R so that 

M> sup ||z(i)|| 

t€[-h,0] 

R = sup{ |/(t,a;)| | t G [-/i,0], ||a;|| < M }. 

For each n choose a continuous function x n : [—h, 0] — > R™ such that = for all i ^ i?„ and 

sup ||a;„(i)|| < Af 
te[-h,o] 

This is possible, for example, by linear interpolation. Now we have, for the continuous function x n 

J f(t,x n (t))dt = K + J (f(t,x n (t)) - f(t,z(t)f)dt 

= K + J (f(t,x n (t))-f(t,z(t)))dt 
< K + 4Ra(k- l)/n 

This proves the desired result. □ 

The following lemma states that when the argmin z f(t, z) is piecewise continuous in t, then we have the desired 
result. 

Lemma 4. Suppose f : [—h,0] x K™ — > R and f/ie hypotheses of Lemma 3 hold. Then the following are equivalent, 
(i) For allye C[-h,0] 

/ f{t,y(t))dt>0. 

J-h 
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(ii) There exists g : [—h,0] — > M which is piecewise continuous and satisfies 



f(t, z) + g(t) > for all t, z 

f g{t)dt = 0. 
J-h 



Proof. Again we only need to show that (i) implies (ii). Suppose (i) holds, then 

,o 

inf / f(t,y(t))dt> 



and hence by Lemma 3 we have 
where r : [— h, 0] — > M™ is given by 



)dt > 



/ r(t), 
-(i) = inf/(i,ar) for all i. 



The function r is continuous on H c since / is continuous on H c x K". Hence by Lemma 2, there exists g such that 
condition (ii) holds, as desired. □ 

We now specialize the result of Lemma 4 to the case of quadratic functions. It is shown that in this case, that under 
certain conditions, the argmin z f(t, z) is piecewise continuous. 

Theorem 5. Suppose M : [—h, 0] — ► is piecewise continuous, and there exists e > such that for all t £ [—h, 0] 

we have 

M 22 (t) > el 



where M is partitioned as 

Mn M12 
_M 21 M 22 

with M 22 : [—h,0] — ► S™. Then the following are equivalent. 



M = 



(i) For all x E R m and continuous y : [—h, 0] — ► R" 





X 


T 


X 




f 




M(t) 




dt>0 


l-h 









(6) 



(ii) There exists a function T : [— h, 0] — > § m which is piecewise continuous and satisfies 



M(t) + 



T(t) 




f 

J-h 



> for all t e [-h, 0] 
T{t) dt = 



Proof. Again we only need to show (i) implies (ii). Suppose x G K", and define 

-. T 



f(t,z) = 



M(t) 



for all t, z. 
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Since by the hypothesis M22 has a lower bound, it is invertible for all t and its inverse is piecewise continuous. 
Therefore z(t) — —M22{t)~ 1 M2i(t)x is the unique minimizer of f(t,z) with respect to z. By the hypothesis (i), we 
have that for all y G C[— h, 0] 

f(t,y(t))dt>0. 
Hence by Lemma 4 there exists a function g such that 

g(t) + f(t, z) > for all t,z 
\dt = 0. 



J-h 



(7) 



The proof of Lemma 2 gives one such function as 

ff(t) = -/(«,*(*)) + £ J f{s,z{a))dt. 



We have 



and therefore <?(£) is a quadratic function of x, say = x T T(t)x, and T : [—h, 0] — > S m is continuous on H c . Then 
equation (7) implies 

-, t 

x 



x 1 T(t)x 



M{t) 



> for all t, z, x 



as required. 



□ 



Notice that the strict positivity assumption on M22 in Theorem 5 is implied by the existence of an e > such that 

V(x) > e\\x\\l 

where ||-||2 denotes the L 2 -norm. 

We have now shown that the convex cone of functions M such that the first term of (3) is nonnegative is exactly 
equal to the sum of the cone of pointwise nonnegative functions and the linear space of functions whose integral is 
zero. Note that in (6) the vectors x and y are allowed to vary independently, whereas (3) requires that x = y(0). It 
is however straightforward to show that this additional constraint does not change the result, using the technique 
in the proof of Lemma 3. 

The key benefit of this is that it is easy to parametrize the latter class of functions, and in particular when M is a 
polynomial these constraints are semidefinite representable constraints on the coefficients of M. 



4 Lie Derivatives 

In this section we will take the opportunity to define the relationship between the functions M and N which define 
the Lyapunov function V, and the functions D and E, which define the Lie derivative of the Lyapunov function V. 



4-1 Single Delay Case 

We first present the single delay case, as it will illustrate the formulation in the more complicated case of several 
delays. Suppose that V E X is given by (3), where M : [-h, 0] ->■ § 2n and N : [-h, 0] x [-h, 0] -> R nxn . Since there 
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is only one delay, if the system is exponentially stable then there always exists a Lyapunov function of this form 
with continuous functions M and TV. Then the Lie derivative of V is 



V(4>) 





T 




<t>{-h) 


D{s) 


4>(-h) 


J>(s) 




<j>(s) 



ds 



hJ-h 



(s) T E(s,t)(j>{t) dsdt. 



Partition D and M as 



M(t) = 



Mn M 12 (i) 
M 21 {t) M 22 {t) 



D(t) 



D u D 12 (t) 
D 21 (t) D 22 (t) 



so that Mu G S™ and Du G S 2 ™. Without loss of generality we have assumed that Mn and Dn are constant. The 
functions D and E are linearly related to M and AT by 



Dn = 



D 12 (t) 



E(s,t) 



A$M n + M n A MnAi 

AjMn 

M 12 (0)+M 21 (0) -M 12 (-/i) 

-M 21 (-/i) 

M 22 (0) 
-M 22 (-ft) 

^M 12 (t)-M 12 (t)+7V(0,t) 

.4f M 12 (t) - JV(-M) 

#22 (t) - -M 22 (i) 



dN(s,t) dN(s,t) 



ds 



dt 



4-2 Multiple- delay case 

We now define the class of functions under consideration for the Lyapunov functions. Define the intervals 

f[-fti,0] if * = 1 

\[-hi,-hi-i) if i = 2, . . . , k. 

For the Lyapunov function V, define the sets of functions 

Yi = [ M : [-ft, 0] -> § 2n I 
Mn (t) = Mn (a) 
M is C 1 on ffj 

F 2 = { TV : [-ft, 0] x [-ft, 0] - 
N(s,t) =N(t,s) T 
N is C 1 on ^ x Hj 



for all s, t e [-ft, 0] 
for alH = 1, . . . , fc I 

for all s, t <G [-ft, 0] 
for alii, j = 1, . . . , k 
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and for its derivative, define 



Z x = {D : [-ft,0] -» §( fc+2 )« | 

Dy(t) = £>«(*) 

D is C° on ffj 
Z 2 = { E : [-ft, 0] x [-ft, 0] -» S" 

£(M) = #(i,s) T 
£ is C° on a x Hi 



for all s, t e [-ft, 0] 
for z, j = 1,.. .,3 

for all i = 1 , . . . , k j 



for all s, i e [-ft, 0] 
for all i,j = 1, . . . , k j 



Here M G Yi is partitioned according to 



M(t) 



Mu Mi 2 (t) 
M 2 i(t) M 22 (t) 



(8) 



where M u e S™ and # € Zi is partitioned according to 



D(t) = 



#n 


#12 


#13 


D 14 (t) 


#21 


#22 


#23 


#24 (t) 


#31 


#32 


#33 


#34 (t) 



#41 (t) #42 (*) #43(0 #44 (t) 



(9) 



where #11, #33, #44 G S n and L> 22 G S^ 1 )™. Let Y = Y 1 x Y 2 and Z = Z 1 x Z 2 . Notice that if M e Fi, then 
M need not be continuous at hi for 1 < i < k — 1, however, we require it be right continuous at these points. We 
also define the derivative M(t) at these points to be the right-hand derivative of M. We define the continuity and 
derivatives of functions in Y 2 , Z\ and Z 2 similarly. 



We define the jump values of M and N at the discontinuities as follows. 



AM(hi) = lim M(t) - lim M(t) 

t— >(— t—*(—hi) — 



for each i = 1, . . . ,k — 1, and similarly define 



AN(h h t)= lim N(s,t) - lim iV(s,t) 
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Definition 6. Define the map L : Y -» Z by (D, E) = L(M, N) if for all t,se [-h, 0] we have 

D n = Aq Mn + Mn A, 



and 



D 12 



+ ^(Mi 2 (0) + M 2 i(0) + M 22 (0)) 



Mix Ax ■ ■ • M u An 

AMi 2 (/ii) ■■• AM 12 (h k -i] 



1 



1 



MiiA fe -Mi 2 (-/i)) 



D 22 
£ 23 = 

#33 = 



- diag(-AM 22 (/ii), . . . , -AM 22 (/i fe _i)) 



--M 22 (-/i) 



Di 4 (t) - iV(0, t) + AjAfiaCi) - Mi 2 (f) 
AiV(-/n,t)+AfMx 2 (t) 

£>24(*) = ; 

AN{-h k ^i,t) + Al_ 1 M X2 {t) 
D ai (t) = AlM 12 (t)-N(-h,t) 
D 44 (t) = -M 22 {t) 



E(s,t) = 



dN(s,t) BN(s,t) 



ds dt ' 

Here M is partitioned as in (8), D is partitioned as in (9), and the remaining entries are defined by symmetry. 

The map L is the Lie derivative operator applied to the set of functions specified by (3); this is stated precisely 
below. Notice that this implies that L is a linear map. 

Lemma 7. Suppose M £ Y\ and N £ Y 2 and V is given by (3). Let (D,E) = L(M,N). Then the Lie derivative of 
V is given by 



V{4>) 



-h 



'<t>(-h )~ 


T 


'<K-ho) 




D(s) 




<t>(-h k ) 




<j>{-hk) 









ds 



(f>(s) T E(s,t)<p(t) dsdt. 



(10) 



-h J-h 



Proof. The proof is straightforward by differentiation and integration by parts of (3). □ 
5 Polynomial Matrices 

In this paper we use piecewise polynomial matrices as a conveniently parametrized class of functions to represent 
the functions M and TV defining the Lyapunov function (3) and its derivative. Theorem 5 has reduced nonnegativity 
of the first term of (3) to pointwisc nonnegativity of a piecewise polynomial matrix in one variable. 

We first make some definitions which we will use in this paper; some details on polynomial matrices may be found 
in Scherer and Hoi [l3| and Kojima Q. We consider polynomials in n variables. As is standard, for a £ N n define 
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the monomial in n variables x a by x a — x" 1 • • • x™ n . We say M is a real polynomial matrix in n variables if for some 
finite set W c N" we have 

M(x) = ^ A a x a 

where A a is a real matrix for each a G W. A convenient representation of polynomial matrices is as a quadratic 
function of monomials. Suppose z is a vector of monomials in the variables x, such as 



z(x) 



1 

2 

•2- 1 <2<2 



For convenience, assume the length of z is d + 1. Let Q G S™( d + 1 ) be a symmetric matrix. Then the function M 
defined by 

M(x) = (/„ ® z(x)) T Q(/„ ® z(x)) (11) 

is an n x n symmetric polynomial matrix, and every real symmetric polynomial matrix may be represented in this 
way for some monomial vector z. If we partition Q as 



Q 



Qi 



Qi 



Qnl ■ ■ ■ Q 



where each Q i} G R( d + 1 ) x ( d + 1 ) ; then the i,j entry of M is 



Mij(x) = z{x) T Q^z{x). 

Given a polynomial matrix M, it is called a sum of squares if there exists a vector of monomials z and a positive 
semidefinite matrix Q such that equation (11) holds. In this case, 

M (x) t for all x 

and therefore the existence of such a Q is a sufficient condition for the polynomial M to be globally pointwise positive 
semidefinite. A matrix polynomial in one variable is pointwise nonnegative semidefinite on the real line if and only 
if it is a sum of squares; see Choi, Lam, and Reznick (lj. Given a matrix polynomial M(x), we can test whether it 
is a sum-of-squares by testing whether there is a matrix Q such that 



M{x) = (J„ <g) z(x)) T Q{I n ® z{x)) 



(12) 



where z is the vector of all monomials with degree half the degree of M. Equation (12) is interpreted as equality of 
polynomials, and equating their coefficients gives a finite set of linear constraints on the matrix Q. Therefore to find 
such a Q we need to find a positive semidefinite matrix sub ject to linear constraints, and this is therefore testable 
via semidefinite programming. See Vandenberghe and Boyd [16[ for background on semidefinite programming. 

5.1 Piecewise polynomial matrices 

A matrix- valued function M : [— h, 0] — > S n is called a piecewise polynomial matrix if for each i = l,...,k the 
function M restricted to the interval Hi is a polynomial matrix. We represent such piecewise polynomial matrices 
as follows. Define the vector of indicator functions g : [—ft, 0] — > R k by 



9i{t) = 



1 if t G Hi 
otherwise 
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for alH = 1, . . . , k and all t G [—ft, 0]. Let z(t) be the vector of monomials 



z(t) 



and for convenience also define the function Z n ^ d : [-ft, 0] — > K nfc ( d + 1 ) x ™ by 

= g(t) ® I n ® z(t). 

Then it is straightforward to show that M is a piecewise matrix polynomial if and only if there exist matrices 
Qi G §™( d+1 ) for i = 1, . . . , k such that 

M(t) =Z n , d (i) T diag(Qi,...,Q fc )Z n , d (t). (13) 

The function M is pointwise positive semidefinite, i.e., 

M{t) y for all t 6 [-ft, 0] 

if there exists positive semidefinite matrices Qi satisfying (13). We refer to such functions as piecewise sum of squares 
matrices, and define the set of such functions 

X n ,d = { Zl d {t)QZ n4 {t) I 

Q = diag(Q 1 ,...,Q fe ), Q t e^ d+1 \Q t hO}- 

If we are given a function M : [—ft, 0] — > S" which is piecewise polynomial and want to know whether it is piecewise 
sum of squares, then this is computationally checkable using semidefinite programming. Naturally, the number of 
variables involved in this task scales as kn 2 (al + l) 2 when the degree of M is 2d. 



5.2 Piecewise Polynomial Kernels 



We consider functions N of two variables s, t which we will use as a kernel in the quadratic form 

/ (j)(s) T N(s,t)(i){t)dsdt (14) 

-h J-h 

which appears in the Lyapunov function (3). A polynomial in two variables is referred to as a binary polynomial. 
A function N : [—ft, 0] x [—ft, 0] — ► S" is called a binary piecewise polynomial matrix if for each i, j G {1, . . . , k} the 
function N restricted to the set Hi x Hj is a binary polynomial matrix. It is straightforward to show that N is a 
symmetric binary piecewise polynomial matrix if and only if there exists a matrix Q E §" fc ( d + 1 ) suc h that 

N(s,t) = Zl d (s)QZ n , d (t). 

Here d is the degree of N and recall 

Z n ,d(t)=g(t)®I n ®z(t). 

We now proceed to characterize the binary piecewise polynomial matrices N for which the quadratic form (14) is 
nonnegative for all <p G C([— ft, 0], R"). We first state the following Lemma. 
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Lemma 8. Suppose z is the vector of monomials 



z(t) 



1 

t 

t 2 



and the linear map A : C[0, 1] — * R d+1 is given by 



A(f>= z{t)4>{t)dt 
Jo 



Then rank A = d + 1 . 



Proof. Suppose for the sake of a contradiction that rank^L < d + 1. Then ranged is a strict subset of R d+1 and 
hence there exists a nonzero vector q 6 R d+1 such that q _L range A. This means 



q T z(t)cj)(t) dt = 



for all (f) e C[0, 1]. Since q T z and </> are continuous functions, define the function v : [0, 1] — ► M by 



u(f) = / q T z(s) ds for all t £ [0, 1]. 
Jo 

Since v is absolutely continuous, we have for every <fi £ C[0, 1] that 

l ,i 

4>{t)dv{t) = I q T z{t)cj>{t) dt 
o Jo 
= 

where the integral on the left-hand-side of the above equation is the Stieltjes integral. The function v is also of 
bounded variation, since its derivative is bounded. The Riesz representation theorem [121 ] implies that if v is of 
bounded variation and 

4>{t)dv{t) = o 

for all <f> € C[0, 1], then v is constant on an everywhere dense subset of (0, 1). Since v is continuous, we have v is 
constant, and therefore q T z(t) = for all t. Since q T z is a polynomial, this contradicts the statement that q ^ 0. □ 



We now state the positivity result. 

Theorem 9. Suppose N is a symmetric binary piecewise polynomial matrix of degree d. Then 



rO 



{s) T N{s,t)(l)(t) dsdt>0 



(15) 



hJ-h 



for all <t> S C([-h, 0],R") if and only if there exists Q E § nk (d+i) suc h t h at 

N(s,t) = Zl d (s)QZ n<d (t) 
QhO, 
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Proof. We only need to show the only if direction. Suppose N is a symmetric binary piecewise polynomial matrix. 
Let d be the degree of N. Then there exists a symmetric matrix Q such that 

N(s,t) = Zl d (s)QZ ntd (t). 

Now suppose that the inequality (15) is satisfied for all continuous functions cj>. We will show that every such Q is 
positive semidefinite. To see this, define the linear map J : C([-h, 0], R") — » R" fe ( d+1 ) by 



Then 



-h J-h 



J<f>= [ {q{t)d)I n ®z{t))<t>{t)dt. 
J-h 



4>{s) T N{s,t)(j){t)dsdt = {J<j)) T Q{J(j)). 



The result we desire holds if rank J = nk(d + 1), since in this case range J = R" fc ( d+1 ). Then if Q has a negative 
eigenvalue with corresponding eigenvector q, there exists <f> such that q = J(f> so that the quadratic form will be 
negative, contradicting the hypothesis. 

To see that rank J = nk(d + 1), define for each % = 1, . . . , k the linear map Li : C[Hi] — > R n by 

L l( j>= j z(t)4>{t)dt 

J Hi 



Then if we choose coordinates for 6 such that 



I -Hi 

\h 2 



□ 



where restriction of <j) to the interval Hj, then we have in these coordinates that J is 

J = diagfXi, . . . , Lk) <8> /„. 
Further, by Lemma 8 the maps Lj each satisfy rankLj = + 1. Therefore rank J = nfc(rf + 1) as desired. 

The following corollary gives a tighter degree bound on the representation of N. 

Corollary 10. Let N be a binary piecewise polynomial matrix of degree 2d which is positive in the sense of Equa- 
tion (15), then there exists a Q E §" fc ( d+1 ) such that 

N(s,t) = Zl d (s)QZ ntd (t) 
QhO, 



Proof. The binary representation used in Theorem 9 had the form 

N(s,t) = Zl 2d (s)PZ nM (t) 
PhO, 



14 



where P G §nfe(2d+i) However, in any such a representation, it is clear that Pij,ij = for i = d + 2, . . . , 2d + 1 
and j = 1, . . . , kn. Therefore, since P >z these rows and columns are and can be removed. Define Q to be the 
reduction of P. Z n ^ is the corresponding reduction of Z n ^d- Then Q G §,nk{d+i) ^ q ^ , and 

N(s,t) = Zl 2d (s)PZ nM (t) = Zl d (s)QZ nid (t). 

□ 



For convenience, we define the set of symmetric binary piecewise polynomial matrices which define positive quadratic 
forms by 

V n4 = {Zl d {s)QZ n4 {t) Q G §" fe(d+1) ,Q b 0}. 



As for £„.(2, if we are given a binary piecewise polynomial matrix N : [— h, 0] x [— h, 0] — > § n of degree 2<i and 
want to know whether it defines a positive quadratic form, then this is computationally checkable using semidefinite 
programming. The number of variables involved in this task scales as (nk) 2 (d + l) 2 . 



6 Stability Conditions 



Theorem 11. Suppose there exist d G N and piecewise matrix polynomials M, T, N, D, U, E such that 



M + 



-D + 



T 



U 





G S 



2n,d 



G S 



(k+2)n,d 



i: 
s: 



N G T n4 
-E G T n ^d 
(D, E) = L(M 7 N) 

T(s)ds = 
U(s) ds = 

Mn y o 

Dn ~< 



Then the system defined by equations (1) and (2) is exponentially stable. 



Proof. Assume M, T, N, D,U, E satisfy the above conditions, and define the function V by (3). Then Lemma 7 
implies that V is given by (10). The function V is the sum of two terms, each of which is nonnegative. The first 
is nonnegative by Theorem 5 and the second is nonnegative since N G T„^. The same is true for V. The strict 
positivity conditions of equations (4) hold since Mn >~ and — D\\ >- 0, and Theorem 1 then implies stability. □ 



The feasibility conditions specified in Theorem 11 are semidcfinitc-representablc. In particular the condition that a 
piecewise polynomial matrix lie in S is a set of linear and positive semidefinitcness constraints on its coefficients. 
Similarly, the condition that T and U integrate to zero is simply a linear equality constraint on its coefficients. Stan- 
dard semidefinite programming codes may therefore be used to efficiently find such piecewise polynomial matrices. 
Most such codes will also return a dual certificate of infcasibility if no such polynomials exist. 
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As in the Lyapunov analysis of nonlinear systems using sum-of-squares polynomials, the set of candidate Lyapunov 
functions is parametrized by the degree d. This allows one to search first over polynomials of low degree, and increase 
the degree if that search fails. 

There are various natural extensions of this result. The first is to the case of uncertain systems, where we would like 
to prove stability for all matrices Ai in some given semialgebraic set. This is possible by extending Theorem 11 to 
allow Lyapunov functions which depend polynomially on unknown parameters. A similar approach may be used to 
check stability for systems with uncertain delays. Additionally, stability of systems with distributed delays defined 
by polynomial kernels can be verified. It is also straightforward to extend the class of Lyapunov functions, since it 
is not necessary that each piece of the piecewise sums-of-squares functions be nonnegative on the whole real line. 
To do this, one can use techniques for parameterizing polynomials nonnegative on an interval; for example, every 
polynomial p(x) — /(x) — (x — l)(x — 2)g(x) where / and g are sums of squares is nonnegative on the interval [1, 2]. 



7 Numerical Examples 



In this section we present the results of some example computations using the approach described above. The 
computations were performed using Matlab, together with the SOSTOOLS [Hi toolbox and SeDuMi pj] code for 
solving semidefinite programming problems. 



7. 1 Illustration. 



Consider the process of proving stability using the results of this paper. The following system has well-known stability 
properties. 

x(t) = -x(t-l) (16) 
A Matlab implementation of the algorithm in this paper has been developed and is available online. This implemen- 
tation returns the following Lyapunov function for system (16). For symmetric matrices, sub-diagonal elements are 
suppressed. 



V(x) = 
where 

and 

Positivity is proven using the function 
and the sum of squares functions 

and 
where 

and 



i: 


x(0)~ 


T 

M(s) 


x(0)~ 


ds + 


LI 




x(s)_ 


x(s)_ 







x(s) T R(s,t)x(t) dsdt 



M(s) = 



27.3 -16.8 + 2.74s 
24.3 + 8.53s 



R(s,t) = 9.08. 
<(s) = -.915+ 1.83s 
Q(s) 

V(s) = Z(s) T LZ(s), 
Z(s) 



13 -3.3 
12.2 



> 



1 s 
1s 



L = 



28.215 5.585 -16.8 -1.973 
13 1.413 -3.3 
24.3 10.365 
12.2 



> 0. 
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This is because — s(s + 1) > for s G [— 1, 0] and 

M(s) 



t(s) 




-s(s + l)Q(s) + V(s) 



Furthermore 

R(s) = 9.08 > 0. 

Therefore, by Theorems 11 and 5, the Lyapunov function is positive. 
The derivative of the function is given by 



V(x) = |° 



" x(0) " 


T 


' x(0) ' 




x(-l) 


D(s) 


x(-l) 


ds 


x(s) _ 









where 

9.3 7.76 -6.34 
-D(s) = 15.77 -7.72 + 2.74s 

8.53 

Negativity of the function is proven using the function 



U(s) 



.0055 + . Oils -.272 -.544s 
-.458 - .916s 



where 



and the sum-of-squares functions 



U(s)ds = 0, 



X(s) 



8.86 1.90 -3.23 
11.54 -3.71 
7.74 



and 
where 



Y(s) = Z{s) T LZ(s), 



Z(s) = 



1 s 
1 s 
Is 



and 



9.3 4.43 7.76 .61896 -6.34 -1.0371 
8.86 1.281 1.9 -2.1929 -3.23 
15.77 5.77 -7.72 .01171 
11.54 -.98171 -3.71 
8.53 3.87 
7.74 



> 0. 
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Negativity follows since 

'U(s) O" 


Therefore, by Theorem 11, the derivative of the Lyapunov function is negative. Stability follows by Theorem 1. 



-D(s) 



-s(s + l)X(s)+Y{s). 



7.2 A single delay. 

We consider the following instance of a system with a single delay. 



x(t) = 



1 
-2 0.1 



x(t) 





1 



x(t - h) 



For a given h, we use semidcfinitc programming to search for a Lyapunov function of degree d that proves stability. 
Using a bisection search over h, we determine the maximum and minimum h for which the system may be shown to 
be stable. These are shown below. 



d 


^rnin 


^max 


1 


.10017 


1.6249 


2 


.10017 


1.7172 


3 


.10017 


1.71785 



When the degree d = 3, the bounds h m i n and /i max are tight Q. For comparison, we include here the bounds obtained 
by Gu using a piecewise linear Lyapunov function with n segments. 



n 




^max 


1 


.1006 


1.4272 


2 


.1003 


1.6921 


3 


.1003 


1.7161 



7.3 Multiple delays. 

Consider the system with two delays below. 



x(t) 



1 
-1 ± 



x(t) 




-1 



x(t-§) + 





1 



x(t - h) 



As above, using a bisection search over h, we prove stability for the range of h below. 



d 


^min 


^max 


1 


.20247 


1.354 


2 


.20247 


1.3722 



Here for degree d = 2, the bounds obtained are tight. Again we include here bounds obtained by Gu Q using a 
piecewise linear Lyapunov function with n segments. 



n 




^max 


1 


.204 


1.35 


2 


.203 


1.372 
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8 Summary 



In this paper we developed an approach for computing Lyapunov functions for linear systems with delay. In general 
this is a difficult computational problem, and certain specific classes of this problem are known to be NP-hard [l5| . 
However, the set of Lyapunov functions is convex, and this enables us to effectively test nonemptiness of a subset of 
this set. Specifically, we parameterize a convex set of positive quadratic functions using the set of polynomials as a 
basis, and the main results here are Theorems 5 and 9. Combining these results with the well-known approach using 
sum-of-squares polynomials allows one to use standard semidefinite programming software to compute Lyapunov 
functions. This gives a nested sequence of computable sufficient conditions for stability of linear delay systems, 
indexed by the degree of the polynomial. In principle this enables searching over increasing degrees to find a Lyapunov 
function, although further work is needed to enhance existing semidefinite programming codes to make this more 
efficient in practice. 

It is possible that Theorem 5 and its proof techniques are applicable more widely, specifically to stability analysis 
of nonlinear systems with delay, as well as to controller synthesis for such systems. One specific extension that is 
possible is analysis of delay systems with uncertain parameters, for which sufficient conditions for existence of a 
Lyapunov function may be given using convex relaxations. It is also possible to analyze stability of nonlinear delay 
systems, in the case that the dynamics are defined by polynomial delay-differential equations. Further extensions to 
allow synthesis of stabilizing controllers are of interest, and may be possible. 
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